# Replication R script for HATRICC QCA
# Created by: Meghan Lane-Fall
# Last updated: 1/14/2023

# Corresponding author: Meghan Lane-Fall
# Corr. auth e-mail: laneme@pennmedicine.upenn.edu
# Corr. auth ORCID: 0000-0001-7050-0017

# GETTING STARTED #################################################
rm(list = ls()) # Clear environment
cat("\014") # Clear console

# STEP 1: INPUT & READ DATA, RENAME VARIABLES#################

# Set working directory
setwd('C:/Users/lanef/Box Sync/Active projects/HATRICC QCA/Meghan R Files/For Harvard Dataverse')

# Load necessary packages
library(QCA)
library(SetMethods)

# Import QCA data for HATRICC
hatQCA <- read.csv("HATRICC_QCAdataset_replication.csv", header = TRUE, row.names=2)


# Rename variables to make the names shorter
colnames(hatQCA)[which(names(hatQCA) == "NEW2ICU")] <- "N2I.u"
colnames(hatQCA)[which(names(hatQCA) == "ICU1")] <- "ICU.u"
colnames(hatQCA)[which(names(hatQCA) == "FIDELITY")] <- "F.u"
colnames(hatQCA)[which(names(hatQCA) == "ANESPRESENT")] <- "AP.u"
colnames(hatQCA)[which(names(hatQCA) == "SURGPRESENT")] <- "SP.u"
colnames(hatQCA)[which(names(hatQCA) == "RNPRESENT")] <- "RP.u"
colnames(hatQCA)[which(names(hatQCA) == "OPPRESENT")] <- "OP.u"
colnames(hatQCA)[which(names(hatQCA) == "ACTORSUM")] <- "SUM.u"
colnames(hatQCA)[which(names(hatQCA) == "TCOMP")] <- "TC.u"
colnames(hatQCA)[which(names(hatQCA) == "FAMILHIGH")] <- "FAMH.u"
colnames(hatQCA)[which(names(hatQCA) == "QUIET")] <- "Q.u"
colnames(hatQCA)[which(names(hatQCA) == "ATTNHIGH")] <- "AH.u"
colnames(hatQCA)[which(names(hatQCA) == "ENGAGEHIGH")] <- "EH.u"
colnames(hatQCA)[which(names(hatQCA) == "PROFHIGH")] <- "PH.u"
colnames(hatQCA)[which(names(hatQCA) == "TEAMHIGH")] <- "TH.u"

#Creating secondary outcome (low fidelity)
hatQCA$FL.u <- 10-hatQCA$F.u

# Outcome names key (*.u is uncalibrated, *.f is fuzzy, *.c is crisp)
# F.u: Fidelity (uncalibrated)
####### FH.f: Fidelity high, fuzzy-set
# FL.u: Fidelity low (uncalibrated)
####### FL.f: Fidelity low, fuzzy-set

# Condition names key
# N2I.u: New to ICU (i.e., not re-admitted; uncalibrated)
####### N2I.c: New to ICU, crisp-set

# AP.u: Anesthesiologist/anesthetist present (uncalibrated)
####### (did not calibrate this condition, unused)

# SP.u: Surgeon present (uncalibrated)
####### (did not calibrate this condition, unused)

# RP.u: RN (nurse) present (uncalibrated)
####### (did not calibrate this condition, unused)

# OP.u: Ordering provider present (uncalibrated)
####### OP.c: Ordering provider present, crisp-set

# SUM.u: Sum of clinicians present, out of 4 possible (uncalibrated)
####### (did not calibrate this condition, unused)

# TC.u: Team complete (i.e., 4/4 clinicians present; uncalibrated)
####### (did not calibrate this condition, unused)

# FAMH.u: Familiarity high (uncalibrated)
####### FAMH.f: Familiarity high, fuzzy-set

# Q.u: Quiet context (uncalibrated)
####### Q.c: Quiet context, crisp-set

# AH.u: Attention high (uncalibrated)
####### AH.f: Attention high, fuzzy-set

# EH.u: Engagement high (uncalibrated)
####### EH.f: Engagement high, fuzzy-set

# PH.u: Professionalism high (uncalibrated)
####### PH.f: Professionalism high, fuzzy-set

# TH.u: Teamwork high (uncalibrated)
####### TH.f: Teamwork high, fuzzy-set


# STEP 2: CALIBRATION #####################################

# Creating empty data frames for calibrated variables (conditions and outcomes)
# 60 empty rows created because I have 60 observations (after deleting obs with missing data)
hatQCA1 <- data.frame(matrix(ncol = 0, nrow = 60)) # for primary outcome
hatQCA2 <- data.frame(matrix(ncol = 0, nrow = 60)) # for secondary outcome

# Calibrating outcomes -----------------------------------------
#Calibrating the primary outcome (high fidelity) as fuzzy-set with direct calibration
hatQCA1$FH.f <- calibrate(hatQCA$F.u, type="fuzzy", method="direct", c(5.5, 6.5, 8.5))
hatQCA1$FH.f <- round(hatQCA1$FH.f, digits = 2) # rounding to 2 places

#Calibrating secondary outcome (low fidelity) as fuzzy-set with direct calibration
hatQCA2$FL.f <- calibrate(hatQCA$FL.u, type="fuzzy", method="direct", c(2.5, 3.5, 4.5))
hatQCA2$FL.f <- round(hatQCA2$FL.f, digits = 2) # rounding to 2 places

# Calibrating conditions (including candidate conditions) ----------------------

#Calibrating condition 1 of 9: N2I.c (crisp-set)
hatQCA1$N2I.c <- calibrate(hatQCA$N2I.u, type = "crisp", method = "direct", c(0.5))
hatQCA2$N2I.c <- calibrate(hatQCA$N2I.u, type = "crisp", method = "direct", c(0.5))

#Calibrating condition 2 of 9: ICU.c (crisp-set)
hatQCA1$ICU.c <- calibrate(hatQCA$ICU.u, type = "crisp", method = "direct", c(0.5))
hatQCA2$ICU.c <- calibrate(hatQCA$ICU.u, type = "crisp", method = "direct", c(0.5))

#Calibrating condition 3 of 9: OP.c (crisp-set)
hatQCA1$OP.c <- calibrate(hatQCA$OP.u, type = "crisp", method = "direct", c(0.5))
hatQCA2$OP.c <- calibrate(hatQCA$OP.u, type = "crisp", method = "direct", c(0.5))

#Calibrating condition 4 of 9: FAMH.f (fuzzy-set with manual calibration)
hatQCA1$FAMH.f <- NA # creates an empty vector
hatQCA1$FAMH.f[hatQCA$FAMH.u>=3] <- 1 # upper threshold
hatQCA1$FAMH.f[hatQCA$FAMH.u<3 & hatQCA$FAMH.u>1] <- 0.6 # more in than out
hatQCA1$FAMH.f[hatQCA$FAMH.u<=1] <- 0 # lower threshold

hatQCA2$FAMH.f <- NA # creates an empty vector
hatQCA2$FAMH.f[hatQCA$FAMH.u>=3] <- 1 # upper threshold
hatQCA2$FAMH.f[hatQCA$FAMH.u<3 & hatQCA$FAMH.u>1] <- 0.6 # more in than out
hatQCA2$FAMH.f[hatQCA$FAMH.u<=1] <- 0 # lower threshold

#Calibrating condition 5 of 9: AH.f (fuzzy-set with manual calibration)
hatQCA1$AH.f <- NA # creates an empty vector
hatQCA1$AH.f[hatQCA$AH.u>=3] <- 1 # upper threshold
hatQCA1$AH.f[hatQCA$AH.u<3 & hatQCA$AH.u>1] <- 0.6 # more in than out
hatQCA1$AH.f[hatQCA$AH.u<=1] <- 0 # lower threshold

hatQCA2$AH.f <- NA # creates an empty vector
hatQCA2$AH.f[hatQCA$AH.u>=3] <- 1 # upper threshold
hatQCA2$AH.f[hatQCA$AH.u<3 & hatQCA$AH.u>1] <- 0.6 # more in than out
hatQCA2$AH.f[hatQCA$AH.u<=1] <- 0 # lower threshold

#Calibrating condition 6 of 9: EH.f (fuzzy-set with manual calibration)
hatQCA1$EH.f <- NA # creates an empty vector
hatQCA1$EH.f[hatQCA$EH.u>=3] <- 1 # upper threshold
hatQCA1$EH.f[hatQCA$EH.u<3 & hatQCA$EH.u>1] <- 0.6 # more in than out
hatQCA1$EH.f[hatQCA$EH.u<=1] <- 0 # lower threshold

hatQCA2$EH.f <- NA # creates an empty vector
hatQCA2$EH.f[hatQCA$EH.u>=3] <- 1 # upper threshold
hatQCA2$EH.f[hatQCA$EH.u<3 & hatQCA$EH.u>1] <- 0.6 # more in than out
hatQCA2$EH.f[hatQCA$EH.u<=1] <- 0 # lower threshold

#Calibrating condition 7 of 9: TH.f (fuzzy-set with manual calibration)
hatQCA1$TH.f <- NA # creates an empty vector
hatQCA1$TH.f[hatQCA$TH.u>=3] <- 1 # upper threshold
hatQCA1$TH.f[hatQCA$TH.u<3 & hatQCA$TH.u>1] <- 0.6 # more in than out
hatQCA1$TH.f[hatQCA$TH.u<=1] <- 0 # lower threshold

hatQCA2$TH.f <- NA # creates an empty vector
hatQCA2$TH.f[hatQCA$TH.u>=3] <- 1 # upper threshold
hatQCA2$TH.f[hatQCA$TH.u<3 & hatQCA$TH.u>1] <- 0.6 # more in than out
hatQCA2$TH.f[hatQCA$TH.u<=1] <- 0 # lower threshold

#Calibrating condition 8 of 9: PH.f (fuzzy-set with manual calibration)
hatQCA1$PH.f <- NA # creates an empty vector
hatQCA1$PH.f[hatQCA$PH.u>=3] <- 1 # upper threshold
hatQCA1$PH.f[hatQCA$PH.u<3 & hatQCA$PH.u>1] <- 0.6 # more in than out
hatQCA1$PH.f[hatQCA$PH.u<=1] <- 0 # lower threshold

hatQCA2$PH.f <- NA # creates an empty vector
hatQCA2$PH.f[hatQCA$PH.u>=3] <- 1 # upper threshold
hatQCA2$PH.f[hatQCA$PH.u<3 & hatQCA$PH.u>1] <- 0.6 # more in than out
hatQCA2$PH.f[hatQCA$PH.u<=1] <- 0 # lower threshold

#Calibrating condition 9 of 9: Q.c (crisp-set)
hatQCA1$Q.c <- calibrate(hatQCA$Q.u, type = "crisp", method = "direct", c(0.5))
hatQCA2$Q.c <- calibrate(hatQCA$Q.u, type = "crisp", method = "direct", c(0.5))

# ANALYSIS #1 - Primary outcome as fuzzy set (FH.f)

# STEP 3 (hatQCA1): NECESSITY #########################################

QCAfit(hatQCA1[, 2:10], hatQCA1$FH.f, necessity = TRUE) # outcome
QCAfit(hatQCA1[, 2:10], 1- hatQCA1$FH.f, necessity = TRUE) # NON-outcome

# STEP 4 (hatQCA1): TRUTH TABLE CONSTRUCTION ##########################

# Create a truth table with final set of conditions:
TThatQCA1 <- truthTable(hatQCA1, "FH.f", 
                        conditions = c("N2I.c", 
                                       "OP.c",
                                       "AH.f",
                                       "Q.c"),
                        complete = TRUE, show.cases = TRUE, 
                        incl.cut = 0.75, sort.by = "incl, n")
TThatQCA1
TThatQCA1$tt # creating object with just minimized truth table
write.csv(TThatQCA1$tt, "HatQCA Truth Table 1.csv") # writing to csv

# Minimizing, conservative solution
sol.hatQCA1cons <- minimize(TThatQCA1, include = "1", details = TRUE, 
                                use.tilde = TRUE)
sol.hatQCA1cons

# Minimizing, parsimonious solution
sol.hatQCA1pars <- minimize(TThatQCA1, include = "?", details = TRUE, 
                                use.tilde = TRUE)
sol.hatQCA1pars

# Examine the simplifying assumptions that were used in the parsimonious solution
sol.hatQCA1pars$SA

# Intermediate Solution
sol.hatQCA1int <- minimize(TThatQCA1, outcome = "FH.f", include = "1, ?", 
                           row.dom = TRUE, details = TRUE, show.cases = TRUE, 
                           use.tilde = TRUE, exclude = c(10), 
                           method = "CCubes")
sol.hatQCA1int


# ANALYSIS #2 - Secondary outcome: low fidelity (FL.c)

# STEP 3 (hatQCA2): NECESSITY ######################################

QCAfit(hatQCA2[, 2:10], hatQCA2$FL.f, necessity = TRUE) # outcome, secondary
QCAfit(hatQCA2[, 2:10], 1- hatQCA2$FL.f, necessity = TRUE) # NON-outcome, secondary

# STEP 4 (hatQCA2): TRUTH TABLE CONSTRUCTION ##########################

TThatQCA2 <- truthTable(hatQCA2, "FL.f", 
                        conditions = c("N2I.c",
                                       "OP.c", 
                                       "EH.f",
                                       "TH.f",
                                       "Q.c"),
                        complete = TRUE, show.cases = TRUE, incl.cut = 0.75, 
                        sort.by = "incl, n")
TThatQCA2
TThatQCA2$tt # creating object with just minimized truth table
write.csv(TThatQCA2$tt, "HatQCA Truth Table 2.csv") # writing to csv

# Minimizing, conservative solution
sol.hatQCA2cons <- minimize(TThatQCA2, include = "1", details = TRUE, 
                            use.tilde = TRUE)
sol.hatQCA2cons

# Minimizing, parsimonious solution
sol.hatQCA2pars <- minimize(TThatQCA2, include = "?", details = TRUE, 
                            use.tilde = TRUE)
sol.hatQCA2pars

# Examine the simplifying assumptions that were used in the parsimonious solution
sol.hatQCA2pars$SA

# Robustness tests

rob.calibrange(raw.data = hatQCA, calib.data = hatQCA1,
               test.cond.raw = "F.u", test.cond.calib = "FH.f",
               test.thresholds = c(5.5, 6.5, 8.5), step = 0.5, max.runs = 1000,
               outcome  = "FH.f", conditions = c("N2I.c", "OP.c", "AH.f", "Q.c"),
               incl.cut = 0.75, n.cut = 1, include = "?")


# Show which cases are in which solution path
settable <- compute("OP.c*AH.f + N2I.c*OP.c*Q.c + N2I.c*AH.f*Q.c", 
                                    data=hatQCA1, separate = TRUE)
settable$handoffnum <- hatQCA$Handoffnum # adding case IDs to table
settable$fidelity <- hatQCA1$FH.f # adding fidelity to table
write.csv(settable, "settable.csv") # writing to csv

# Figure 1: Outcome

# Figure 1a: raw fidelity data
par(mfrow = c(2, 2), lwd = 2)

ggplot(hatQCA, aes(x=F.u, xmin = 0)) +
  geom_bar(fill = "Dark Gray", col = "Steel Blue", lwd = 1) +
  labs(x='Fidelity, # steps of 10 (outcome)', 
       y='Frequency') +
  scale_x_continuous(limits = c(0, 10.5), breaks = c(0, 2, 4, 6, 8, 10))

# Figure 1b: raw vs. calibrated fidelity data
calibgraph <- data.frame(raw = hatQCA$F.u, fuzzy = hatQCA1$FH.f)

ggplot(calibgraph, aes(x=raw, y=fuzzy)) +
  geom_point(shape = 21, 
             fill = "Dark Gray",
             col = "Steel Blue",
             stroke = 1.5, 
             size = 4) +
  labs(x='Fidelity (raw score)', 
       y='Fidelity (fuzzy score)')

# Figure 2: conditions

# Figure 2a: New to ICU
CatVar <- within(hatQCA, {
  N2I.cat <- NA
  N2I.cat[hatQCA$N2I.u < 0.5] <- "No (ICU re-admission)"
  N2I.cat[hatQCA$N2I.u > 0.5] <- "Yes"
})

ggplot(CatVar, aes(x=N2I.cat)) + 
  geom_bar(fill = "Dark Gray", col = "Steel Blue", lwd = 1) +
  labs(x='New to ICU', y='Number of cases with condition')

# Figure 2b: Ordering provider present
CatVar <- within(hatQCA, {
  OP.cat <- NA
  OP.cat[hatQCA$OP.u < 0.5] <- "Absent"
  OP.cat[hatQCA$OP.u > 0.5] <- "Present"
})

ggplot(CatVar, aes(x=OP.cat)) + 
  geom_bar(fill = "Dark Gray", col = "Steel Blue", lwd = 1) +
  labs(x='Ordering provider present', y='Number of cases with condition')

# Figure 2c: Attention high
CatVar <- within(hatQCA, {
  AH.cat <- NA
  AH.cat[hatQCA$AH.u < 2] <- "Unsatisfactory"
  AH.cat[hatQCA$AH.u > 1 & hatQCA$AH.u < 3] <- "Satisfactory"
  AH.cat[hatQCA$AH.u >2] <- "Superior"
})

CatVar$AH.cat <- factor(CatVar$AH.cat, levels = 
                          c("Unsatisfactory", "Satisfactory", "Superior"))

ggplot(CatVar, aes(x=AH.cat)) + 
  geom_bar(fill = "Dark Gray", col = "Steel Blue", lwd = 1) +
  labs(x='Attention rating', y='Number of cases with condition')

# Figure 2d: Quiet environment

CatVar <- within(hatQCA, {
  Q.cat <- NA
  Q.cat[hatQCA$Q.u < 0.5] <- "Absent"
  Q.cat[hatQCA$Q.u > 0.5] <- "Present"
}) 

ggplot(CatVar, aes(x=Q.cat)) + 
  geom_bar(fill = "Dark Gray", col = "Steel Blue", lwd = 1) +
  labs(x = 'Quiet environment', y='Number of cases with condition')

############ END #################